Introduction to Portfolio Optimization

Piotr Lipiński
Computational Intelligence Research Group, Institute of Computer Science, University of Wroclaw, Poland
lipinski@cs.uni.wroc.pl

Abstract:

This notebook illustrates the concept of portfolio optimization in the context of the financial time series from the Warsaw Stock Exchange.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

import datetime
import zipfile

%matplotlib inline

# import plotly.offline as py
# import plotly.figure_factory as ff
# import plotly.graph_objs as go

# py.init_notebook_mode(connected=True)

# py_config = {'displayModeBar': False, 'showLink': False, 'editable': False}

# from ipywidgets import interact, interactive, fixed, interact_manual, widgets

STOCK_QUOTATIONS_ARCHIVE_FILE_NAME = 'mstall.zip'
STOCK_NAMES_FILE_NAME = 'WIG20.txt'

def load_stock_quotations(stock_names, filename):
    s = {}
    with zipfile.ZipFile(filename) as z:
        for stock_name in stock_names:
            with z.open(stock_name + '.mst') as f:
                s[stock_name] = pd.read_csv(f, index_col='<DTYYYYMMDD>', parse_dates=True)[['<OPEN>', '<HIGH>', '<LOW>', '<CLOSE>', '<VOL>']]
                s[stock_name].index.rename('time', inplace=True)
                s[stock_name].rename(columns={'<OPEN>':'open', '<HIGH>':'high', '<LOW>':'low', '<CLOSE>':'close', '<VOL>':'volume'}, inplace=True)
    return pd.concat(s.values(), keys=s.keys(), axis=1)

1. Notowania spółek giełdowych

stock_names = pd.read_csv(STOCK_NAMES_FILE_NAME, index_col=0, names=['stock name']).index.values
stock_names.sort()
number_of_stocks = len(stock_names)

stock_quotations = load_stock_quotations(stock_names, STOCK_QUOTATIONS_ARCHIVE_FILE_NAME)
stock_quotations.fillna(method='ffill', inplace=True)
stock_quotations.tail()

	ALIOR					ASSECOPOL					...	SYNTHOS					TAURONPE
	open	high	low	close	volume	open	high	low	close	volume	...	open	high	low	close	volume	open	high	low	close	volume
time
2018-05-08	67.95	68.05	65.70	66.05	420855.0	44.60	44.84	44.22	44.48	57018.0	...	4.9	4.91	4.88	4.89	267394.0	2.28	2.29	2.23	2.25	6635586.0
2018-05-09	66.10	68.55	65.90	68.30	420927.0	44.30	45.00	44.30	45.00	47211.0	...	4.9	4.91	4.88	4.89	267394.0	2.26	2.35	2.25	2.35	5186793.0
2018-05-10	69.90	73.50	69.20	72.55	521518.0	44.98	45.00	44.68	45.00	31980.0	...	4.9	4.91	4.88	4.89	267394.0	2.34	2.44	2.34	2.44	6516527.0
2018-05-11	72.30	73.45	72.00	72.80	360528.0	44.98	45.28	44.78	44.84	106236.0	...	4.9	4.91	4.88	4.89	267394.0	2.44	2.44	2.38	2.44	3747375.0
2018-05-14	72.80	73.05	72.35	72.80	156323.0	45.16	45.30	44.80	44.94	33078.0	...	4.9	4.91	4.88	4.89	267394.0	2.44	2.45	2.35	2.35	4641902.0

5 rows × 100 columns

stock_name = 'TAURONPE'

plt.figure(figsize=(12, 4))
plt.plot(stock_quotations[stock_name]['close'].dropna())
plt.xlabel('time')
plt.ylabel('price')
plt.title(stock_name)
plt.show()

2. Stopy zwrotu spółek giełdowych

Rozpatrujemy $d$ spółek giełdowych, $\mathcal{A}_1, \mathcal{A}_2, \ldots, \mathcal{A}_d$. Niech $v^i_t$ oznacza wartość akcji spółki $\mathcal{A}_i$ w dniu $t$ (dla celów obliczeniowych, przyjmujemy, że wartość akcji jest określana przez cenę zamknięcia). Stopą zwrotu spółki $\mathcal{A}_i$ w dniu $t$ nazywamy $$r^i_t = \frac{v^i_t - v^i_{t-1}}{v^i_{t-1}}.$$

Niech $t_0$ będzie ustalonym przyszłym dniem. Zakładamy w związku z tym, że nie mamy informacji o notowaniach spółek giełdowych w dniu $t_0$, ale znamy notowania z dni poprzednich. Stopę zwrotu spółki $\mathcal{A}_i$ w dniu $t_0$ traktujemy więc jako zmienną losową, oznaczmy ją przez $R_i$, i możemy określić oczekiwaną stopę zwrotu $\mathbb{E}[R_i]$ oraz wariancję stopy zwrotu $\mathbb{Var}[R_i]$ estymując je na okresie ostatnich $\Delta t$ dni, czyli na okresie $[t_0 - \Delta t, t_0 - 1]$. Analogicznie, oznaczając przez $\mathbf{R} = (R_1, R_2, \ldots, R_d)^T$ wektor losowy stóp zwrotu wszystkich rozpatrywanych spółek, możemy określić kowariancje $\mathbb{Cov}[\mathbf{R}]$ i korelacje $\mathbb{Corr}[\mathbf{R}]$ stóp zwrotu. Zatem:

• oczekiwaną stopę zwrotu $\mathbb{E}[R_i]$ spółki $\mathcal{A}_i$ w dniu $t_0$ definiujemy jako średnią stopą zwrotu z ostatnich $\Delta t$ dni, czyli z okresu $[t_0 - \Delta t, t_0 - 1]$,

• wariancję $\mathbb{Var}[R]$ spółki $\mathcal{A}_i$ w dniu $t_0$ definiujemy jako wariancję na ostatnich $\Delta t$ dniach, czyli na okresie $[t_0 - \Delta t, t_0 - 1]$,

• zależności między spółkami w dniu $t_0$ określamy przez kowariancję/korelację ich stóp zwrotu na ostatnich $\Delta t$ dniach, czyli na okresie $[t_0 - \Delta t, t_0 - 1]$.

---

UWAGA: Przez dzień rozumiem tutaj dzień roboczy, a nie dzień kalendarzowy. Prowadzi to do pewnych nieścisłości (stopa zwrotu dla poniedziałku to de facto trzydniowa stopa zwrotu), bo tak liczone stopy zwrotu różnią się od stóp zwrotu liczonych zazwyczaj dla instrumentów pozbawionych ryzyka, m.in. lokat bankowych, gdzie oprocentowanie liczy się dla wszystkich dni kalendarzowych.

---

Poniższe przykładowe obliczenia ilustrują estymowanie oczekiwanej stopy zwrotu i wariancji dla ostatniego dnia zarejestrowanych notowań.

delta_t = 90
stock_returns = stock_quotations.xs('close', level=1, axis=1).pct_change()
stock_returns_m = stock_returns[-delta_t-1:-1].mean()
stock_returns_s = stock_returns[-delta_t-1:-1].std()
stock_covariances = stock_returns[-delta_t-1:-1].cov()
stock_correlations = stock_returns[-delta_t-1:-1].corr()

2.1. Oczekiwana stopa zwrotu i ryzyko (określone przez wariancję stopy zwrotu) rozpatrywanych spółek

plt.figure(figsize=(12, 4))
plt.bar(np.arange(stock_returns_m.size)+0.1, stock_returns_m, 0.8)
plt.gca().xaxis.set_tick_params(width=0)
plt.xticks(np.arange(stock_returns_m.index.size)+0.5, stock_returns_m.index.values, rotation='vertical')
plt.xlim([0, stock_returns_m.size])
plt.title('expected return rate')
plt.show()

plt.figure(figsize=(12, 4))
plt.bar(np.arange(stock_returns_s.size)+0.1, stock_returns_s**2, 0.8)
plt.gca().xaxis.set_tick_params(width=0)
plt.xticks(np.arange(stock_returns_s.index.size)+0.5, stock_returns_s.index.values, rotation='vertical')
plt.xlim([0, stock_returns_s.size])
plt.title('variance of return rate')
plt.show()

2.2. Korelacja stóp zwrotu rozpatrywanych spółek

plt.figure(figsize=(9, 6))
plt.imshow(stock_correlations, cmap='hot', interpolation='none')
plt.colorbar()
plt.gca().xaxis.set_tick_params(width=0)
plt.xticks(np.arange(stock_correlations.columns.size), stock_correlations.columns.values, rotation='vertical')
plt.gca().yaxis.set_tick_params(width=0)
plt.yticks(np.arange(stock_correlations.columns.size), stock_correlations.columns.values, rotation='horizontal')
plt.title('correlation of return rates')
plt.show()

2.3. Spółki dominujące, zdominowane i niezdominowane

Zastanawiając się w którą spółkę giełdową zainwestować kapitał, rozważamy dwa parametry: oczekiwaną stopę zwrotu inwestycji oraz ryzyko inwestycji (definiowane przez wariancję stopy zwrotu). Rozpatrując zbiór spółek giełdowych, można wprowadzić pojęcia spółki dominującej i spółki zdominowanej w następujący sposób:

• Mówimy, że spółka giełdowa $\mathcal{A}$ dominuje spółkę giełdową $\mathcal{B}$, jeżeli oczekiwana stopa zwrotu spółki $\mathcal{A}$ jest nie niższa niż spółki $\mathcal{B}$ oraz ryzyko spółki $\mathcal{A}$ jest niższe niż ryzyko spółki $\mathcal{B}$ .
• Spółka $\mathcal{B}$ jest spółką zdominowaną, jeśli istnieje spółka $\mathcal{A}$, która ją dominuje.
• Spółka $\mathcal{B}$ jest spółką niezdominowaną, jeśli nie istnieje spółka $\mathcal{A}$, która ją dominuje.

is_dominated = np.zeros(stock_returns_m.size)
for i in range(stock_returns_m.size):
    for j in range(stock_returns_m.size):
        if (i != j) and (stock_returns_m[i] <= stock_returns_m[j]) and (stock_returns_s[i] > stock_returns_s[j]):
            is_dominated[i] = 1
            break

plt.figure(figsize=(9, 6))
plt.scatter(stock_returns_s, stock_returns_m, color='#000000')
plt.xlabel('standard deviation of return rate')
plt.ylabel('expected return rate')
for i in range(stock_returns_m.size):
    if is_dominated[i] == 0:
        name = stock_returns_m.index.values[i]
        plt.annotate(
            name,
            xy = (stock_returns_s[name], stock_returns_m[name]), xytext = (40, 40),
            textcoords = 'offset points', ha = 'right', va = 'bottom', size = 8,
            bbox = dict(boxstyle = 'round, pad=0.5', fc = 'yellow', alpha = 0.5),
            arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3, rad=0'))
plt.show()

3. Teoria portfela Markowitza

Zamiast zastanawiać się w którą spółkę giełdową zainwestować cały kapitał, można zastanawiać się jak podzielić rozpatrywany kapitał na części i zainwestować je w poszczególne spółki (stworzyć portfel inwestycji).

3.1. Wprowadzenie

Niech $\mathbf{R} = (R_1, R_2, \ldots, R_d)^T$ oznacza wektor losowy stóp zwrotu wszystkich rozpatrywanych spółek.

Niech $\boldsymbol{\mu} = (\mu_1, \mu_2, \ldots, \mu_d)^T = \mathbb{E}[\mathbf{R}] \in \mathbb{R}^d$ oznacza wektor oczekiwanych stóp zwrotu rozpatrywanych spółek.

Niech $\boldsymbol{\Sigma} = \mathbb{Cov}[\mathbf{R}] \in \mathbb{R}^{d \times d}$ oznacza macierz kowariancji stóp zwrotu rozpatrywanych spółek.

Zakładamy, że rozpatrywane spółki są instrumentami obarczonymi ryzykiem, czyli wariancja ich stóp zwrotu jest niezerowa (dokładniej: zakładamy, że macierz kowariancji $\boldsymbol{\Sigma}$ jest dodatnio określona - skąd wynika, że jest macierzą odwracalną).

Portfelem nazywamy wektor $\mathbf{p} = (p_1, p_2, \ldots, p_d)^T \in \mathbb{R}^d$, taki że $\sum_{i=1}^d p_i = 1$, który definiuje rozbicie kapitału na części do zainwestowania w poszczególne spółki ($p_i$ to wielkość kapitału do zainwestowania w spółkę $\mathcal{A}_i$).

Stopa zwrotu portfela $\mathbf{p}$ jest więc zmienną losową $R_\mathbf{p} = \mathbf{p}^T \mathbf{R}$.

Oczekiwana stopa zwrotu portfela wynosi zatem $\mathbb{E}[R_\mathbf{p}] = \mathbb{E}[\mathbf{p}^T \mathbf{R}] = \mathbf{p}^T \mathbb{E}[\mathbf{R}] = \mathbf{p}^T \boldsymbol{\mu}$.

Ryzyko portfela, określone przez wariancję jego stopy zwrotu, wynosi zatem $\mathbb{Var}[R_\mathbf{p}] = \mathbb{Var}[\mathbf{p}^T \mathbf{R}] = \mathbf{p}^T \mathbb{Cov}[\mathbf{R}] \mathbf{p} = \mathbf{p}^T \boldsymbol{\Sigma} \mathbf{p}$.

UWAGA: Powyższe podejście teoretyczne nie pokrywa się z podejściem obliczeniowym (wartość oczekiwana i wariancja stopy zwrotu portfela estymowana na podstawie jego wartości w ustalonym okresie nie pokrywają się z wartościami liczonymi według powyższych wzorów). Pozornie, dla ustalonego dnia $t_0$, oczekiwaną stopę zwrotu portfela (i analogicznie jej wariancję), można byłoby określić estymując je na okresie ostatnich $\Delta t$ dni, czyli na okresie $[t_0 - \Delta t + 1, t_0]$. Portfel definuje jaka część całkowitego kapitału (domyślnie równego 1) jest inwestowana w poszczególne spółki: w $i$-tą spółkę inwestowane jest $p_i$, czyli kupowane jest $\frac{p_i}{v^i_{t_0}}$ akcji tej spółki. Zatem, całkowita wartość takiego portfela w dniu $t$ wynosiła $$\sum_{i=1}^d \frac{p_i}{v^i_{t_0}} v^i_t,$$ czyli stopa zwrotu takiego portfela w dniu $t$ wynosiła $$\frac{\sum_{i=1}^d \frac{p_i}{v^i_{t_0}} v^i_t - \sum_{i=1}^d \frac{p_i}{v^i_{t_0}} v^i_{t-1}}{\sum_{i=1}^d \frac{p_i}{v^i_{t_0}} v^i_{t-1}} = \frac{\sum_{i=1}^d \frac{p_i}{v^i_{t_0}} (v^i_t - v^i_{t-1})}{\sum_{i=1}^d \frac{p_i}{v^i_{t_0}} v^i_{t-1}} = \frac{\sum_{i=1}^d p_i r^i_t \frac{v^i_{t-1}}{v^i_{t_0}}}{\sum_{i=1}^d p_i \frac{v^i_{t-1}}{v^i_{t_0}}} \approx \frac{\sum_{i=1}^d p_i r^i_t}{\sum_{i=1}^d p_i} = \sum_{i=1}^d p_i r^i_t,$$ o ile $\frac{v^i_{t-1}}{v^i_{t_0}}$ jest bliskie $1$. Jeśli $t$ byłoby równe $t_0 + 1$, to oczywiście równość zachodziłaby. W przeciwnym przypadku możemy otrzymać rozbieżność między bezpośrednio estymowaną oczekiwaną stopą zwrotu portfela, a oczekiwaną stopą zwrotu portfela wyliczoną z modelu zawierającego estymowane statystyki stóp zwrotu rozpatrywanych spółek. Rozbieżność ta wynika z faktu, że próbujemy przekształcić portfel wyrażający procentowy udział poszczególnych spółek w całkowitym kapitale portfela na liczbę akcji poszczególnych spółek, a ze względu na zmiany cen akcji w czasie, ustalenie konkretnej liczby akcji w wybranym dniu $t_0$ może powodować zmiany procentowego udziału poszczególnych spółek w całkowitym kapitale portfela w czasie. Innymi słowy, w rozpatrywanej wyżej estymacji, w dniach $t \neq t_0$ analizowaliśmy portfel nie będący portfelem $\mathbf{p}$!

3.2. Optymalizacja portfela

Portfel optymalny to portfel minimalizujący ryzyko dla ustalonej oczekiwanej stopy zwrotu. Można go formalnie określić jako rozwiązanie problemu optymalizacji funkcji $$\frac{1}{2} \mathbf{p}^T \boldsymbol{\Sigma} \mathbf{p}$$ z ograniczeniami $$\mathbf{p}^T \boldsymbol{\mu} = e, \quad \mathbf{p}^T \mathbf{1} = 1,$$ gdzie $e$ jest ustaloną oczekiwaną stopą zwrotu portfela.

Stosując metodę mnożników Lagrange'a, można rozwiązać ten problem analitycznie otrzymując rozwiązanie [1]: $$\mathbf{p} = \frac{1}{D} ((e C - A) \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu} + (B - e A) \boldsymbol{\Sigma}^{-1} \mathbf{1}) = \frac{1}{D} (B \boldsymbol{\Sigma}^{-1} \mathbf{1} - A \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}) + e \frac{1}{D} ((C \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu} - A \boldsymbol{\Sigma}^{-1} \mathbf{1}),$$ gdzie $$A = \boldsymbol{\mu}^T \boldsymbol{\Sigma}^{-1} \mathbf{1}, \qquad B = \boldsymbol{\mu}^T \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}, \qquad C = \mathbf{1}^T \boldsymbol{\Sigma}^{-1} \mathbf{1}, \qquad D = B C - A^2.$$

Rozwiązania problemu optymalizacji portfela można przedstawić jako $$\mathbf{p} = \mathbf{p}_1 + e \mathbf{p}_2,$$ gdzie $$\mathbf{p}_1 = \frac{1}{D} (B \boldsymbol{\Sigma}^{-1} \mathbf{1} - A \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}), \qquad \mathbf{p}_2 = \frac{1}{D} ((C \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu} - A \boldsymbol{\Sigma}^{-1} \mathbf{1}).$$

d = stock_returns_m.size

Sinv = np.linalg.inv(stock_covariances)

A = stock_returns_m.T.dot(Sinv.dot(np.ones(d)))
B = stock_returns_m.T.dot(Sinv.dot(stock_returns_m))
C = np.ones(d).T.dot(Sinv.dot(np.ones(d)))
D = B * C - A**2

p1 = 1/D * (B * Sinv.dot(np.ones(d)) - A * Sinv.dot(stock_returns_m))
p2 = 1/D * (C * Sinv.dot(stock_returns_m) - A * Sinv.dot(np.ones(d)))

p = np.array([p1 + 0.0001 * i * p2 for i in range(1000)]).T

p_m = p.T.dot(stock_returns_m)
p_s = np.sqrt(np.diag(p.T.dot(stock_covariances.dot(p))))

plt.figure(figsize=(18, 12))

plt.subplot(2, 2, 1)
plt.scatter(stock_returns_s**2, stock_returns_m, color='#000000')
plt.plot(p_s**2, p_m, '-r')
plt.xlabel('variance of return rate')
plt.ylabel('expected return rate')

plt.subplot(2, 2, 2)
plt.scatter(stock_returns_s**2, stock_returns_m, color='#000000')
plt.plot(p_s**2, p_m, '-r')
plt.xlim([0, 0.0016])
plt.ylim([-0.01, 0.01])
plt.xlabel('variance of return rate')
plt.ylabel('expected return rate')

plt.subplot(2, 2, 3)
plt.scatter(stock_returns_s, stock_returns_m, color='#000000')
plt.plot(p_s, p_m, '-r')
plt.xlabel('standard deviation of return rate')
plt.ylabel('expected return rate')

plt.subplot(2, 2, 4)
plt.scatter(stock_returns_s, stock_returns_m, color='#000000')
plt.plot(p_s, p_m, '-r')
plt.xlim([0, 0.04])
plt.ylim([-0.01, 0.01])
plt.xlabel('standard deviation of return rate')
plt.ylabel('expected return rate')

plt.show()

4. Teoria portfela w praktyce

W ostatnich latach powstaje wiele rozszerzeń klasycznej teorii portfela, które często prowadzą do problemów optymalizacji wymagających złożonych algorytmów rozwiązywania. Wiele modeli wprowadza inne niż wariancja miary ryzyka (m.in. semiwariancja). Wiele modeli odrzuca też założenie o braku ograniczeń krótkiej sprzedaży (ang. short sale), czyli możliwości posiadania ujemnej liczby akcji (w klasycznej teorii portfela dopuszcza się, że współrzędne $p_i$ portfela mogą być ujemne).

References

[1] Markowitz, H., M., Portfolio Selection, Journal of Finance, vol.7, no.1, 1952, pp.77–91.